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abstract! 
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A  procedure  for  performing  sensitivity  analysis  of  data  vis-a-vis 
the  multivariate  error  model,  with  or  without  structural  variables,  is  pro¬ 
posed.  It  is  based  on  an  estimation  procedure,  a  generalization  of  maximum 
likelihood,  which  is  indexsd  on  a  user-specified  value  c.  The  sensitivity 
analysis  proceeds  in  the  following  way:  the  parameters  of  a  set  of  data 
presumed  to  follow  a  multivariate  normal  error  model  are  estimated  vis  maximum) 
likelihood  (c=0) .  The  parameters  are  next  estimated  for  an  index  cj>0 
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20.  Abstract  (cont'd) 

__ _ The  sensitivity  in  parameter  estimates  is  noted.  The  parameters  are  next 

estimated  for  a  user- specified  index  C2>C}>jtf  and  the  sensitivity  of  the 
parameter  estimates  to  the  change  in  c  is  again  noted.  This  process  may 
have  one  or  more  values  of  c>0.  If  the  parameter  estimates  are  sensitive 
functions  of  o,  the  model  and  the  data  are  not  mutually  consistent  and 
both  require  further  detailed  study.  The  items  of  data  which  are  most 
contributory  to  this  sensitivity  are  identified  by  an  examination  of  obser¬ 
vational  weights  which  are  a  byproduct  of  the  analysis.  If  the  parameters 
and  observational  weights  are  non-sensitive  to  changes  in  the  index  c, 
then  one  can  generally  be  confident  of  the  interned,  consistency  of  the 
data  and  the  error  model.  Fixed  values  of  the  index  c  provide  robust 
estimation  procedures  for  model  parameters.  Asymptotic  relative  efficiency 
and  influence  functions  are  provided  for  fixed  values  of  c.  The  results  of 
a  small  Monte  Carlo  study  suggest  that  the  asymptotic  properties  of  the 
estimators  are  rapidly  attained.  Several  illustrations  are  given. 
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ABSTRACT 


A  procedure  for  performing  sensitivity  analyses  of  data  vis-a-vis 
the  multivariate  error  model,  with  or  without  structural  variables,  is 
proposed.  It  is  based  on  an  estimation  procedure,  a  generalization  of 
maximum  likelihood,  which  is  indexed  on  a  user-specified  value  c.  The 
sensitivity  analysis  proceeds  in  the  following  way:  the  parameters  of 
a  set  of  data  presumed  to  follow  a  multivariate  normal  error  model  are 
estimated  via  maximum  likelihood  (c**0) .  The  parameters  are  next  esti¬ 
mated  for  an  index  c^  >  0.  The  sensitivity  in  parameter  estimates  is 
noted.  The  parameters  are  next  estimated  for  a  user-specified  index 
=2  >  Cl  >  0  and  th.  sensitivity  of  the  pare-ter  estimates  to  the 
change  in  c  is  again  noted.  This  process  may  have  one  or  more  values 
of  c  >  0.  If  the  parameter  estimates  are  sensitive  functions  of  c,  the 
model  and  the  data  are  not  mutually  consistent  and  both  require  further 
de ted. led  study.  The  items  of  data  which  are  most  contributory  to  this 
sensitivity  are  identified  by  an  examination  of  observational  weights 
which  eure  a  byproduct  of  the  analysis .  If  the  parameters  and  observa¬ 
tional  weights  are  non-sensitive  to  changes  in  the  index  c,  then  one  can 
generally  be  confident  of  the  internal  consistency  of  the  data  and  the 
error  model.  Fixed  values  of  the  index  c  provide  robust  estimation 
procedures  for  model  parameters.  Asymptotic  relative  efficiencies  and 
influence  functions  are  provided  for  fixed  values  of  c.  The  results  of 
a  small  Monte  carlo  study  suggest  that  the  asymptotic  properties  of  the 
estimators  are  rapidly  attained.  Several  illustrations  are  given. 
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1.  Introduction 

A  number  of  robust  procedures  for  the  analysis  of  p-dimensional 
Gaussian  (normal )  data  have  been  proposed  by  Gnanadesikan  and 
Kettenring  (1972) ,  Devlin*  Gnanadesikan,  and  Kettenring  (1975) ,  Maronna 
(1976),  Huber  (1977),  Hosteller  and  Tukey  (1977),  Devlin,  Gnanadesikan, 
and  Kettenring  (1981) .  A  number  of  numerical  and  statistical  character¬ 
istics  of  these  procedures  have  been  evaluated  in  Devlin  et  al.  (1981) . 
Summaries  and  further  discussion  of  multivariate  robust  procedures  are 
given  in  Gnanadesikan  (1977) ,  Barnett  and  Lewis  (1978) ,  and  Huber  (1981) . 
Despite  the  impressive  range  of  work  in  the  problem  of  robustness  and 
allied  topics,  there  is  still  a  need  for  a  unified  approach  to  the 
subject,  preferably  one  based  on  an  objective  function. 

He  propose  a  procedure  for  the  analysis  of  multivariate  Gaussian 
data  which  is  based  on  a  generalization  of  maximum  likelihood  which  has 
recently  been  considered  in  the  univariate  context  by  Paulson.  Presser, 
and  Nicklin  (1982).  The  estimation  procedure  depends  on  a  user-specified 
index  c  which  may  be  varied  to  determine  (a) ,  the  response  of  the 
parameter  estimates  to  variation  in  this  index,  and  (b),  the  response 
of  observational  weights  associated  with  each  observation  to  variation 
in  this  index.  The  response  surface  generated  permits  a  sensitivity 
analysis  which  is  useful  in  assessing  the  mutual  consistency  of  the  model 
and  the  data  considered  as  a  single  entity.  If  the  parameter  estimates 
*nd  final  observational  weights  are  not  sensitive  to  variation  in  the 
user-specified  index,  then  it  is  highly  likely  that  the  data  and  the 
model  are  internally  consistent.  If  the  parameter  estimates  and  final 
weights  are  sensitive  to  variation  in  the  user-specified  index,  then 
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it  is  likely  that  tha  data  and  the  nodal  are  not  internally  consistent 
for  any  of  a  variety  of  reasons.  Thus  the  procedure  we  are  proposing  can 
be  quite  useful  in  the  identification  of  potential  outliers. 

Our  procedure  produces  robust  estimators  of  the  mean  vector  and  the 
covariance  matrix  by  simply  taking  a  fixed  value  of  the  user-specified 
index  c.  Asymptotic  properties  are  given  for  several  values  of  the 
user-specified  constant  and  for  several  values  of  the  dimension  p.  Tne 
procedure  is  easy  to  use  and  computationally  attractive.  Several  examples 
and  illustrations  are  provided. 


2.  The  Simultaneous  Estimation  Procedure 

He  suppose  that  x, ,  x„,  . ..,  x  constitutes  a  random  sample  from, 

~~x  ~n 

tentatively,  a  multivariate  Gaussian  density 


f(»|y.  B>  -  | 2iti3 | exp{-1s (g-y)T  e'-ks-y)}  (2.1) 

where  jj  is  a  p*l  vector  of  location  parameters,  and  Q  *  (d^)  is  a  p*p 
positive  definite  matrix  of  covariances..  The  estimators  for  y  and  D,  which 
we  shall  term  self-critical ,  are  determined  from  maximization  of 


l  -  A  l }  fC(gllu'B> 
c  c  i-J  lB(U.B.ollc/ll+0) 


(2.2) 


where 


1+e , 


Q(y#B»c)  *  I  t“  "'(s|y»B)<*« 

'Rp 

{(1+C)P  <27T>CP|B)C>"a* 


(2.3) 
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It  is  easily  shown  that 


lim  1  -  l  «  l  log  f(s.  (y,D), 

c-0  ci-l  1 


the  likelihood  function.  Upon  differentiation  of  (2.2)  with  respect  to  u 
and  0,  and  setting  the  resulting  expressions  to  Q/  we  find  that  the 
estimators  for  y  and  Q  are  determined  from  the  simultaneous  zeros  of 


n  3log  f, 

jifi0{(i+c,  , 


3y 


i  fi°  id+c) 


3log  f, 


li2Sfi}„  q  , 
3Q  w 


(2.4) 


(2.5) 


t-i  - 

where  we  have  suppressed  the  arguments  of  the  density  and  Q  and  written 
f^  for  ftx^Jy.D)  and  Q  for  Q(y,D;c).  Equations  (2.4)  and  (2.5)  may  be 
rearranged  so  that  the  estimators  for  y  and  D,  say  y  and  D,  satisfy  the 
implicit  equations 


n 

l 

i-1  i 


ic 


n 


(2.6) 


?  (x.-Wfr--1 


(1+c) 


i-1 


5i-r''5i-y>  vic 


(2.7) 


i-1 


ic 


where 


vt(j  -  exp(—  (x£-y)TD_1 


(2.8) 


4 


The  weighting  in  (2.6)  and  (2.7)  reflects  a  self -weighting,  namely  by  the 
density  itself.  Thus  the  assumed  density  and  the  data  interact  in  pro¬ 
ducing  parameter  estimators,  and  for  this  reason  we  have  called  the 
resulting  estimators  self-critical. 

The  system  (2.6)- (2. 7)  defines  a  family  of  estimators  indexed  on  the 
user-specified  constant  c.  He  will  presently  address  the  question  of 
how  c  should  be  chosen  for  practical  applications.  It  is  recommended 
that  the  estimators  be  computed  numerically  by  a  fixed  point  algo¬ 
rithm  with  initial  trial  solutions  the  usual  maximum  likelihood 
estimators  of  y  and  D.  Since  the  function  £Q  is  concave  (Rao,  1965, 

pp.  447-449) ,  £  will  also  be  concave  for  at  least  a  neighborhood  of  c  *  0 
c 

and  will,  for  this  neighborhood,  always  have  a  unique  solution.  The 
system  (2.6)- (2.6)  need  not  always  have  a  unique  solution  in  general 
but  we  have  not  encountered  any  difficulties  with  the  estimation  procedure 
in  a  wide  variety  of  experiences  with  real  and  simulated  data  so  long 
as  the  trial  solution  for  c  >  0  was  the  maximum  likelihood  (c*0)  estimators 
and  c  was  not  too  far  from  zero.  As  c  moves  further  and  further  away  from 
zero,  the  procedure  may  break  down  numerically  or  provide  multiple  solutions. 

It  is  easily  verified  that  the  expectations  of  the  left-hand  sides  of 
(2.4)  and  (2.5)  are  the  zero  vector  and  matrix  respectively.  In  fact,  the 
procedure  was  developed  to  guarantee  just  this;  see  Paulson  et  al.  (1982)  for 
further  discussion. 


3 .  Some  Examples 


A  few  examples  at  this  point  will  be  useful  in  showing  the  potential 
of  the  family  of  estimators  (2.6)- (2.8)  in  the  analysis  of  multivariate  data. 
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Example  3.1 

Consider  the  bivariate  data  of  Figure  1  (Anderson  (1958) ) .  The 
central  cloud  of  points  represents  the  original  data  to  which  a  single 
outlier  has  been  added  and  rotated  through  360°  in  9°  increments  starting 
at  the  point  labeled  1.  Table  1  presents  the  results  of  the  tradi¬ 
tional  maximum  likelihood  analysis  applied  to  the  original  data  plus 
the  single  outlier  contrasted  with  the  self-critical  analysis  with  c  = 

.3.  The  effect  of  the  single  outlier  on  the  estimate  of  the  means 

and  y2  is  slight  for  the  c  =  . 3  case.  The  effect  of  the  single  outlier 

2  2 

on  the  estimates  of  the  variances  a,  and  is  substantial  for  c  «  0 

1  2 

but  is  less  than  23%  for  c  =  .3  regardless  of  where  the  outlier  is 
placed.  The  effect  on  the  estimate  of  correlation  for  c  *  0  is  dramatic, 
the  estimate  varying  from  .80  to  .27.  However,  for  c  =>  .30  the  esti¬ 
mate  of  correlation  varies  only  from  .746  and  .800.  In  this  case, 
graphical  methods  make  it  easy  to  see  what  is  happening,  but  when  the 
dimensional i ty  p  increases,  it  becomes  increasingly  difficult  to  iso¬ 
late  problems,  or  more  importantly,  the  elements  of  the  data  which  just 
might  lead  to  the  most  interesting  structure  in  the  data.  The  automatic 
weighting  of  the  observations  via  the  self-critical  procedure  provides 
a  means  at  higher  dimensions  for  indicating  unusual  observations  or 
structure.  The  weight  v ,  for  the  outlying  observations  is  the  lowest 

j  • 

among  all  the  weights. 

Example  3.2 

The  extent  to  which  the  family  of  estimators  (2. 6) -(2. 8)  can  be 
useful  in  identifying  structure  is  now  considered.  Each  of  the  25  3-vector 
observations  obtained  by  Billmeyer  and  Rich  (1978)  from  color  matching  on 
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a  MacAdam  Colorimeter  are  shown  in  Table  2.  The  three  variables  are 
the  percentages  of  red,  blue,  and  green  in  the  final  adjustment  by  the 
observer  in  his  attempt  to  match  a  given  color.  Analysis  of  similar 
data,  obtained  from  22  sessions  with  a  single  observer,  led  Brown  et  al. 
(1956)  to  conclude  that  the  data  were  trivariate  normal.  [The  data  are 
all  ostensibly  from  a  single  trivariate  population.]  However,  normal 
probability  plots  of  each  variable  would  show  that  observations  16,  19, 
20,  21,  22,  and  25  seem  the  most  unusual  in  one  or  more  dimensions  and 
an  Andrews-type  plot  corroborates  this. 

we  applied  equations  (2. 6) -(2. 8)  to  these  data,  starting  first  at 

c  *  0  and  then  progressively  increasing  c  through  .1,  .3,  to  .5.  Several 

observations  sure  immediately  highlighted  by  a  dramatic  change  in  the 

final  weights  w.  .  Severed,  components  of  the  covariance  matrix  also 
3C 

show  marked  sensitivity  to  the  variation  in  c.  As  c  increases  to  .5, 

12  of  the  weights  are  reduced  to  low  values.  The  results  for  c  *  .5 
are  also  presented  in  Table  2.  The  low  weights  w.  _  =  v.  _/Ev . 

j,*->  D  »  •  ->  D  r  •  5 

in  Table  2  pinpoint  observations  14-25  for  further  study.  This  analysis 
led  to  careful  reexamination  of  the  original  data  cards  which  showed 
that  the  first  13  observations  were  those  of  a  single  observer,  A  say, 
while  the  last  12  were  those  of  12  different  observers. 

Figures  2a,  b,  and  c  show  all  two-dimensional  scatter  plots  of  this 
sample  with  observer  A's  observations  labeled  by  A's  and  all  other  ob¬ 
servations  labeled  by  B's.  Greek  letters  represent  various  combinations 
of  observations  which  are  too  close  to  differentiate  in  a  plot  of  this 
scale.  For  example,  £  represents  2  A  observations  and  1  B  observation. 
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TABLE  2 


Sensitivity  Analysis  of  Trivariate  Data 


c*.125  c« . 3  c*. 5 


Point 

Mo. 

X1 

X2 

X3 

Wj » .3 

A 

Wj,.5 

1 

86.07 

58.27 

30.72 

.051 

.064 

0.102 

2 

83.45 

56.50 

29.45 

.053 

.077 

0.115 

3 

82.07 

55.43 

28.71 

.053 

.075 

0.097 

4 

83.58 

56.32 

28.85 

.049 

.060 

0.043 

5 

82.95 

56.21 

30.12 

.050 

.065 

0.053 

6 

83.62 

56.38 

29.51 

.048 

.059 

0.086 

7 

86.25 

58.39 

30.92 

.050 

.063 

0.092 

8 

84.79 

57.62 

30.15 

.047 

.046 

0.050 

9 

82.87 

55.99 

29.87 

.050 

.065 

0.064 

10 

79.74 

53.89 

27.70 

.053 

.055 

0.036 

11 

81.84 

55.26 

29.10 

.051 

.065 

0.091 

12 

84.10 

57.08 

30.08 

.049 

.061 

0.085 

13 

86.24 

58.34 

30.51 

.051 

.061 

0.087 

14 

80.55 

54.32 

26.23- 

.051 

.057 

0.2 (-4) 

15 

81.08 

54.54 

25.32 

.046 

.036 

0.2 (-8) 

16 

101.2 

70.38 

33.97 

•l(-4) 

.1(-18) 

0.6 (-41) 

17 

84.92 

57.75 

28.08 

.033 

.011 

0.3 (-6) 

18 

83.78 

56.52 

26.49 

.044 

.030 

0.1 (-8) 

19 

60.54 

40.34 

20.06 

.020 

.7 (-8) 

0.8 (-14) 

20 

62.33 

41.38 

19.11 

.024 

.2 (-6) 

0.2 (-13) 

21 

72.09 

48.46 

22.50 

.042 

.004 

0.2 (-8) 

22 

63.38 

58.98 

32.66 

.5  (-7) 

.8(-23) 

0.2 (-39) 

23 

81.45 

54.53 

23.78 

.030 

.007 

0.6 (-18) 

24 

78.89 

53.37 

26.21 

.048 

.039 

0.2 (-3) 

25 

92.50 

61.76 

29.02 

.007  •' 

.3  (-5) 

0.1 (-16) 

Figure  2  (continued) 


(c)  X2  v».  X3 


Scatter  Plots  for  Exaople  3.2 
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He  observe  that  the  observations  o£  A  are  tightly  grouped  while  the 
remainder  are  more  generally  scattered.  Thus  the  self-critical  weights 
were  able  to  separate  the  observations  by  a  single  observer  from  those 
of  other  origins,  a  feat  which  the  eye  would  be  hard-pressed  to  match 
without  the  aid  of  the  separate  labels  as  given  on  the  scatter  plot a. 
Indeed,  if  the  observations  were  not  labeled  according  to  their  genesis, 
one  would  be  tempted  to  pick  out  the  B's  lying  at  the  extremes  as  being 
potential  problem  points  and  consequently  miss  the  relevant  structure 
entirely.  This  ability  to  cluster  is  also  implied  by  the  form  of  the 
influence  functions  which  we  present  in  a  later  section. 


4 .  Some  Asymptotics 


>  A  A 

The  estimators  y  and  V  defined  by  (2.4)- (2.8)  are  M-estimators  for 
location  and  the  covariance  matrix  by  virtue  of  (2.4)  and  (2.5).  Further¬ 
more,  if  the  observations  x  are  transformed  according  to  y,  =  a  +  Ax^ , 
then  u  *  a  +  au  and  6  ■  ADAT.  The  estimator  £  is  positive  definite 

*»y  "n  — y  -V 

with  probability  one  as  soon  as  n  >  p.  However,  some  numerical  diffi¬ 
culties  may  surface  if  sufficiently  many  of  the  final  weights  vic  are 
approximately  zero  so  that  there  are  fewer  than  p  v^  which  are  larger  than 
zero.  In  this  latter  cue  V  wil ■  be  algorithmically  singular. 

A  A 

Standard  arguments  may  be  used  to  show  that  the  estimators  Ta  and  Y 
defined  as  the  consistent  zeros  of  (2.4)  and  (2.5)  are  asymptotically 
Gaussian. 


Define 


i  \  *C<*.,IU'D) 
o|  lQ(«,B>c)]c/(1+c) 


(4.1) 
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31  31 

Then  the  jth  score  function  is  — ggl  and  E  (— *  0.  The  self- 
critical  estimators  (Q,g)  are  asymptotically  Gaussian  with  mean  0 
and  variance-covariance  matrix 


So 


h”1  E  H-1 
2q  ~0~0 


where  the  8 , 3*  element  of  H  and  E  are  respectively 

MJ  ~0 


(4.2) 


(4.3) 


(4.4) 


for  any  j  ■  l,2,...,n  and  8,8'  ranges  over  811,1  dll'  d12' 

•  •  *,<iip'd22'  *  •  -d2p'  *  *  .,dpp*  1116  ®*Pr®ss*-ons;  (4.2)-(4.4)  are  evaluated  at 
the  (unknown)  true  values  of  y  and  D,  y^  and  DQ.  When  c  -  0  in  (4.2)- 

(4.4)  the  egressions  are  evaluated  by  a  limiting  argument  and  all  reduce 
to  the  usual  maximum  likelihood  expressions. 

The  covariance  matrix  CQ  may  be  written  as 
A  0 

So  *  (4.5) 

0  B 


where  A  is  the  pxp  covariance  matrix  for  the  estimator  W  and  B  is  the 

**  «w  — 

>1P (p+l)x»jp(p+l)  covariance  matrix  for  the  estimator  D.  The  off-diagonal 
blocks  of  0  indicate  that  the  estimators  y  and  0  are  asymtotically  inde¬ 
pendent.  The  matrix  A  may  be  easily  evaluated  at  the  multivariate  normal 
distribution  and  is  given  by 

A  -  (1+c)P+2(1+c)‘,(p+2)D  (4.6) 

with  elements  a^  ■  covfli^,  y^) .  The  form  of  B  is  considerably  more 
complex  and  we  evaluate  it  only  for  particular  cases. 
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We  now  give  the  general  form  for  the  matrix  B  at  the  multivariate 
normal  distribution  with  mean  vector  0,  covariance  matrix  I.  The 
particular  variances  and  covariances  will  be  indicated  by  the  subscripts 
used.  We  have 


a32  *  (3c2  +  4c+2) (c+l)^*2 (2c+l)  ^P*4^  for 

jj 

Bg2  -  {c+l)P+4(2c+l)“1*{p+4>  ,  for  j  +  kj 

jh 

m  c2  (c+l)^2  (2c+l)  111  ^P4"4^  f  for  j  ^  k; 

jj  Hck 

Bn  a  *  0  »  for  all  other  cases. 

^j^Jt 


(4.7) 


Let  us  consider  now  the  vector  of  estimators  of  covariances 

A  A  A 

dll'  **12'  d22*  At  tile  bivariate  standard  normal  with  correlation  p, 
the  matrix  B  of  (4.5)  is  given  by  (4.2)  with 


and 


50-K2 


H0  "  K1 


3c  +4c+2 


-1 

2p 

-P2 

2p 

-2-2P2 

2p 

2 

-P 

2p 

-1 

(4.8) 


*’2p(3c  +4c+2) 


2  2  2 
c  +2P  (c+1) 


-2p(3c2+4c+2)  4(c+1)2+4P2(2c2+2c+1)  -2p(3c2+4c+2) 

c2+2p2 (c+1) 2  -2P(3c2+4c+2)  3c2+4c+2 


J 


(4.9) 


where  K. 


•j(27T)  c|d|  1,0  2 (c+1)  2  and  K2  ■  %  (2ir)"'2c|D|"c~2(2c+l)_3. 


IS 


The  general  form  of  (4.8)  and  (4.9)  Is  quite  complicated  and  is  not  given 
here.  It  can  be  approximated  in  practice  from  the  asymptotic  results, 
namely  for  (4.2)-(4.4)  evaluated  at  y,  V. 

Table  3  provides  the  asymptotic  efficiencies  of  the  self-critical 
estimators  (c>0)  relative  to  maximum  likelihood  estimators  (c«0)  for  the 

A  A 

estimators  ,  and  d^  for  selected  values  of  the  constant  c.  For  fixed  c, 
there  is  a  decrease  in  the  asymptotic  relative  efficiency  for  the  estimator 

A 

p  as  the  correlation  P  increases.  The  decrease  becomes  more  pronounced 
for  larger  values  of  c.  We  will  see  in  the  sequel  that  the  relative 
efficiencies  of  the  self-critical  estimators  approach  their  asymptotic 
relative  efficiencies  rapidly.  For  samples  as  small  as  n  ■  20,  the  rela¬ 
tive  efficiencies  are  about  equal  to  those  given  in  Table  3. 

Figures  3a,  b,  c  provide  the  asymptotic  relative  efficiences  of  the 
self-critical  estimators  y^,  d^,  p  as  a  function  of  c  and  dimensionality 
p  when  sampling  is  assumed  to  be  from  the  p-variate  circular  normal  dis¬ 
tribution.  For  fixed  c  we  see  that  the  efficiency  decreases  as  p  increases 
and  that  efficiency  decreases  as  c  increases.  For  a  given  value  of  c, 
variance  and  correlation  estimators  are  less  efficient  than  those  of 
location.  It  should  be  emphasized  that  these  figures  pertain  to  the 
p-dimensional  normal  with  covariance  matrix  I  only.  The  efficiencies  for 
covariance  estimators  will  be  somewhat  lower  in  the  general  case  as  is 
indicated  in  Table  3.  Maronna  (1976)  reports  asymptotic  relative  effi¬ 
ciencies  for  the  multivariate  version  of  Huber's  proposal  2  at  the  unit 
spherical  normal  for  Winsorization  proportions  q  »  .2,  .3,  .5,  1.0. 

It  is  interesting  that,  for  fixed  q,  the  efficiencies  increase  with 
increasing  dimensionality. 


Asymptotic  Relative  Efficiency 


.9 


5.  Influence  Functions 


The  influence  function  is  the  aost  important  determinant  of  qual¬ 
itative  robustness  since  many  other  robustness  characteristics  of  an 
estimator  nay  be  derived  from  it.  The  influence  function  characterizes 
the  (asymptotic)  response  of  an  estimator  to  an  additional  observation 


as  a  function  of  where  the  observation  falls.  The  shape  of  the  influ¬ 


n 


ence  function  is  determined  primarily  by  the  score  function 
for  9  m  jj  or  D.  In  order  to  guarantee  that  single  (or  multiple)  observa¬ 
tions  do  not  exert  an  excessive  influence  on  the  response  of  the  esti¬ 
mator  to  single  observations,  we  require  that  the  influence  f miction  be 
bounded.  It  is  also  advantageous,  in  our  opinion,  if  the  influence 
functions  are  re-descending  to  zero. 

The  influence  function  for  jj  and  6  are  both  bounded  and  re-descending 


to  zero  for  c  >  0.  Let  l  denote  Idle  function  when  argument  x.  of  (4.1) 

cx  -3 


is  replaced  by  the  non-random  vector  x.  The  influence  function  for  y 

mt  *0 

at  the  p-dimensional  normal  distribution  with  parameters  y  and  D  is 
given  by 


IC 


/  v  />VV‘ 

’  *  •\^)  V 


(5.1) 


-  (l+c)>5(p+2>  (x-y)  exp(-  |  (x-y)T  j"l(yji)  , 
for  any  j  -  1,1,..., n.  This  influence  funation  is  illustrated  in 
Figaros  4a  and  b  where  we  exhibit  contours  of  the  first 
of  the  vector  (5.1)  with  e  ■  %  and  c  ■  H  respectively  at 
standard  bivariate  normal  distribution  with  sore  correlation. 

These  figures  and  other  figures  not  givsn  hors  indicate  that  to  choose  a 
of  robustness  we  choose  a  larger  value  of  e.  The  contour 
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for  the  first  component  of  the  self-critical  estimator  y  for  the  standard 
bivariate  normal  with  p  =  .9  and  c  ■  **  is  shown  in  Figure  5.  The  nearly 
elliptical  contours  have  the  line  of  perfect  positive  correlation  as  their 
major  axis.  Moderately  distant  observations  that  are  close  to  this  axis 
have  greater  influence  on  the  location  estimator  than  those  that  are 
closer  but  along  the  minor  axis.  Note  also  that  the  component  x2  of  x 
has  an  influence  on  the  estimator  the  first  component  of  • 

The  influence  function  for  the  self-critical  estimator  6  cannot  be 

written  in  simple  form  in  general.  For  the  p  =  2  case  we  illustrate  it 

a  a  a 

in  terms  of  the  vector  of  estimators  w  =  (v^v^  ,v  )  for  tlie  standard 
bivariate  normal  distribution  with  correlation  p .  It  may  be  shown  that 


1C($,  N2(0,£)  = ( 1+c )  ‘ 


f  (1+Ox^  -  1  \ 
(l+c)x1x2  -  P 


C  T  -1 

exp (-  —  x  P  x) 


U+c)x‘  -  1 

\ 

1  P 

l  p  1  . 

c  *  %  for  correlations  p  =  0  and  .9  in  Figures  6a  and  b.  In  Figures 


where  P  = 


We  present  contours  of  this  function  for  d. .  with 
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7a  and  b  we  show  contours  of  the  influence  function  for  the  self-critical 
estimator  of  P  under  the  same  conditions. 

All  the  influence  functions  are  bounded  for  c  >  0  and  are  characterized  by 
closed  contours.  The  closed  contours  indicate  that  the  self-critical 
procedure  will  be  useful  in  clustering  data.  The  clustering  nature  of 
the  estimation  procedures  explains  why  we  were  able  to  identify  the 
13  observations  by  the  same  individual  in  example  3.2.  Just  as  in  this 
example,  the  procedure  will  generally  focus  on  the  tightest  cluster  of 
observations  as  the  parameter  c  is  increased.  The  estimation  procedure 


<b)  P  •  0.9 

Influence  curve  contour*  for  the  self •critical  variance  estimate, 
at  the  bivariate  standard  noraal  distribution,  c  •  % 


<  <o 
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may  also  be  useful  in  classification  problems,  but  we  do  not  pursue  the 
subject  here. 

Devlin  et  al.  (1975)  present  the  influence  functions  for  the  maxi- 
mum  likelihood  estimator  of  P  and  a  5*  trimmed  estimator.  Neither  of 
these  influence  functions  are  bounded,  and  the  trimmed  estimator  is 
unable  to  limit  the  influence  of  observations  in  certain  regions  of 
the  x  space. 

6.  Monte  Carlo  Study 

It  has  been  long  recognized  that  some  alternative  univariate  esti¬ 
mators  such  as  the  at- trimmed  mean  lose  little  efficiency  for  an  exactly 
normal  parent  but  exhibit  much  better  and  more  stable  performance  under 
departures  from  it  (Huber,  1972,  for  example).  To  evaluate  the  perform¬ 
ance  of  the  estimators  of  section  2  with  c  >  0,  a  limited  Monte  Carlo 
study  was  conducted.  Simulated  bivariate  normal  samples  were  used  to 
gauge  losses  in  efficiency  and  slightly  non-normal  samples  were  used  to 
measure  stability  and  performance.  The  performances  of  the  usual  vector 
of  sample  means  and  covariance  matrix,  the  vector  of  10%  (on  each  side) 
trimned  means  combined  with  a  covariance  matrix  estimator  employing  5% 

(on  each  side)  trimmed  variances  and  the  estimators  £  and  V  for  several 
values  of  c  are  examined  under  several  sampling  situations. 

TO  study  the  performance  of  the  selected  estimators  in  sampling  from 

normal  populations,  100  samples  of  size  20  from  each  of  the  bivariate 

2  2 

normal  populations  with  ■  0,  /o2  “  0.5,  1.0,  3.0  and  p  »  0.0, 

0.5,  0.9,  0.99  were  constructed.  The  same  random  number  seed  was  used  to 
generate  samples  for  all  correlations  and  variance  ratios. 
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Tables  of  the  Monte  Carlo  means  and  variances  of  the  population 
parameters  estimated  are  presented  in  Table  4.  In  this  and  the  other 
tables  of  this  section  MEAN  represents  the  sample  mean  vector  or  the 
usual  sample  variances  and  correlations,  TR.  MEAN  represents  the  10% 
on  each  side)  vector  of  trimmed  means  or  the  5%  trimmed  variances  and 
correlations  and  SC  (4) ,  SC  (7-) ,  SC  (7)  represent  the  self-critical  esti- 
mators  of  section  2  with  parameters  c  -  -j  respectively.  All 

estimators  provide  reasonable  estimates  for  location,  but  SC(~)  and 
TR.  Mean  have  slightly  inferior  values  vis-a-vis  the  other  estimators. 

The  c  *  0  and  c  »■  ~r  cases  have  the  smallest  variances.  Comparisons  of 

O 

the  Monte  Carlo  means  and  variances  for  the  variance  ratio  and  corre¬ 
lation  yield  similar  conclusions.  Note  that  the  efficiencies  computed 
by  dividing  the  c  -  0  variances  by  the  c  >  0  variances  are  generally 

good  approximations  to  the  values  given  in  Table  3.  An  exception 
X  A 

occurs  for  c  *«  —  for  the  estimator  p  of  correlation.  An  analogous  char¬ 
acterization  of  the  estimators  is  drawn  if  mean  squared  errors  are  used 
instead  of  variances,  but  we  do  not  present  the  results  here. 

To  study  stability  of  the  estimators  under  slight  departures  from 
a  Gaussian  parent,  bivariate  normal  samples  were  constructed  and  con¬ 
taminants  of  one  of  two  types  were  added  to  the  basic  normal  variates. 

In  the  first  case,  one  simulated  bivariate  Cauchy  was  added  to  standard 
bivariate  Gaussian  samples  of  size  19.  In  the  second  case,  a  mixture  of  a 
random  number  r  of  bivariate  Gaussian  variates  divided  by  a  uniform  deviate 
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and  n  -  r  bivariate  Gaussian  variate  was  constructed.  The  number  r  was 

generated  from  a  binomial  distribution  with  parameters  n  ■  20  and 

probability  .05.  The  parameters  of  the  n  -  r  bivariate  Gaussian  variates 

2  2 

were  *  0,  /02  -  1  and  3,  and  p  -  0,  .5,  .9,  .99.  Table  5 

depicts  the  Monte  Carlo  means  and  variances  for  each  estimator,  based  on 
25  samples  generated  from  each  population.  The  Monte  Carlo  means  for 
location,  variance,  and  correlation  show  that  the  self-critical  estimators 
and  the  trimmed  estimators  retain  their  integrity  with  respect  to  the 
basic  underlying  assumption  of  Gaussianity.  This  is,  however,  not  as 
important  as  the  fact  that  this  stability  concerning  the  basic  cluster 
allows  potential  contaminants  or  model  departures  of  various  types  to 
be  identified. 

A  range  of  experience  with  both  read,  and  simulated  data  corroborates 
the  above  results.  The  estimators  perform  well  in  practice,  especially 
when  used  in  the  context  of  a  sensitivity  analysis. 


Monte  Carlo  Results  for  Non-Normal  Samples  of  Size  20 
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7.  Multivariate  Regression 


He  now  show  that  structured  multivariate  data,  specifically  a  linear 

model,  may  be  naturally  treated  by  the  self -critical  approach  of  section  2 

Let  y.,  y_,  . . . ,  y  be  Nxl  vectors  representing  N  independent  observations 
—1  ~2  —P 

on  each  of  p  correlated  random  variables.  He  assxane  the  model  repre¬ 


sentation 


y.  -  X0j  +  ,  j  «  1,2,...,  p  (7.1) 

where  X  is  N*q  and  considered  given,  §.  is  q*l  and  unknown,  and  is 
an  error  or  disturbance  vector  which  is  distributed  as  For 

a  single  value  of  j ,  (7.1)  is  a  univariate  regression  model.  Define 

y  **  (y. ,  y,,  ...,  XL)/  o  ■  (u. ,  u_,  ...,  u) ,  b  -  (3  ,  8  ,  ...8) 

(Nxp)  12  ^  (Nxq)  1  ~2  (qxp)  ^ 

and  then  (7.1)  may  be  written  as 

Y  -  XB  +  U  .  (7.2) 

He  wish  to  estimate  D  and  B.  Let  y . ,  x. ,  and  u.  denote  the  row  vectors 
of  Y,  X,  and  u  of  dimensions  Ixp,  lxq,  and  lxp  respectively.  Since 
u^  »  B,  the  density  of  the  is  given  by 

5  f<Hll2»8)  -  I21®!”1*  ®*p(-'*(Xi-5i5,2"1(Xi"si5)T  )  (7*3> 

The  objective  function  to  be  maximized  for  B  and  D  is,  from  section  2, 


1  * 

*C  *  i-l[{Q(B,Dic)>c/(1+c) 


3 


On  differentiating  with  respect  to  B  and  D,  setting  the  resulting  ex¬ 
pressions  to  0,  and  solving  we  find  that  B  and  D  satisfy  the  implicit 
equations 

~  (iii-3-  ~1  -  \  *  ~  /) 

(7.5) 


jJjSi1?!  «w(-  f  <?i-Si5>e'l<jr1-ii?>Tj| 


and 


N 


D  -  (1+c)  £  wt (£-x.B)  (y.-x.B)  , 

i«l  ~  ~A 


(7.6) 


where 


f  0,D) 


i  N 

I  fC(y.-x.Bj  0,D) 
i-1  ~i  ~1~  ~  - 


i  ial » 2 1 ... |N. 


(7.7) 


It  is  clear  that  if  c  *  0,  these  equations  reduce  to  the  usual  nmy-tninm 
likelihood  equations.  For  c  >  0,  the  estimators  are  not  equation-by¬ 
equation  univariate  estimators  since  the  covariance  structure  of  D  is  now 
allowed  to  play  a  role  in  the  estimation  process.  It  is  interesting  that, 
from  (7.5)  to  (7.7), 


N 


r  T«  r  T  T  c 

i|1  *i?i  Hi  "  J^iSi  Xi  ~  *i  S 

H  _  H  _  /  N  m  \  ,  N 

*  ~  l  «i?i \  (  l  Si?i  ?i  )  l  *i*i 

i-1  i-1  \  i-1  1  1  V  i-1  1  i 


N 


N 


(7.8) 
T. 


*i 


-  ( 


In  the  same  manner  we  have 

R 


-  5  • 


(7.9) 


■fe’.v 
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The  equations  (7.8)  and  (7.9)  are  weighted  versions  of  the  usual  ortho¬ 
gonality  relations  associated  with  the  maximum  likelihood  or  least 
squares  procedures. 

This  regression  procedure  with  c  >  0  is  useful  in  providing  pro¬ 
tection  against  difficulties  in  the  y  values,  for  example  y  outliers, 
but  is  not  especially  useful  if  there  are  difficulties  associated  with 
the  factor  space,  for  example,  x  outliers.  This  is  easily  seen  by 
rewriting  (7.5)  in  terms  of  the  score  function  for  B.  This  score 
function  is  bounded  and  re-descending  in  the  residuals  -  x^B  but  is 
unbounded  in  the  x^.  Thus  a  single  bad  x^  value  may  ruin  an  entire 
regression  analysis.  An  extension  of  the  approach  we  provide  here  will 
be  useful  in  dealing  with  this  problem,  but  we  do  not  pursue  it  here. 
Belsley,  Kuh  and  Welsch  (1980)  provide  a  nice  discussion  of  the  uni¬ 
variate  version  of  this  problem. 

For  completeness,  we  provide  an  example  in  which  we  compare  the 
regression  equations  obtained  by  maximum  likelihood  and  those  obtained 
from  (7. 5)- (7. 7).  It  is  found  that  the  solutions  are  sensitive  to  the 
change  in  c  from  0  to  j  and  thuB  there  are  some  potential  difficulties 
with  the  data  or  the  model  or  both.  It  is  our  job,  having  been  fore¬ 
warned,  to  find  out  where  the  difficulties  might  be.  However,  this  is 
not  always  possible. 

Example  7.1:  Data  extracted  from  a  study  of  the  effects  of  a  change 
in  environment  on  blood  pressure  is  included  in  Ryan  et  al.  (1976). 

In  this  study,  anthropologists  measured  the  blood  pressure  and  other 
characteristics  of  Peruvian  Indian  males  over  age  21,  who  had  migrated 
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from  primitive  environments  at  high  altitudes  to  modem  lower  altitude 
areas  of  Peru.  Previous  studies  in  Africa  suggested  that  such  migra¬ 
tions  might  increase  blood  pressure  at  first,  but  that  the  blood 
pressure  would  tend  to  return  to  normal  with  time.  The  39  observa- 
tiono  included  in  this  reference  are  listed  in  Table  6?  the  dependent 
variables  are  systolic  and  diastolic  blood  pressure  and  the  independent 
variables  chosen  from  those  available  are  F,  the  fraction  of  life  in 
the  new  environment,  W,  the  weight  in  kilograms,  and  S,  a  skin-fold 
measure  of  general  obesity. 

We  will  consider  the  linear  model  containing  the  aforementioned 
independent  variables  and  a  constant  term.  The  MLE  for  the  full  model  are 


'  55.51 

33.21 

-28.41 

-9.26 

A 

y  — 

85.62 

28.58  ' 

1.39 

0.74 

»  h  * 

28.58 

108.36 

-0.26 

0.001 

4 

. 

The  order  in  the  matrix  of  parameter  estimates  is  constant,  fraction  of 

life,  weight,  and  obesity  measure  with  the  first  column  representing  the 

systolic  blood  pressure  equation  and  the  second  column  representing  the 

diastolic  blood  pressure  equation.  The  self-critical  estimates  with 
1 

parameter  c  -  —  are 
4 


58.59 

34.91 

21.31 

-1.08 

1.28 

0.58 

-0.24 

0.22 

76.98  17.26 

17.26  65.57 

; 


Table  6 

Characteristics  of  Peruvian  Indians  in  Example  7.1 


34 


Heo\40ueH('iegi4NonHNHH4nee<fN(iNNNntennHOMnat'OH 


MHH4nee4NnNNsntennHOMA««H 


09  ■» 

ik 

o  . 


o^Nno'0<®'0HOiflHn>jo'0i>O'0l''NC(io»iflNiA9i^rNaMvinH-jo^ 

r*4  rH  fH  rM  H  H 


SnnNoovO'OinmnaO'Oioifl-jic.foiAMincioONONrsHffi-ioOiOioHWOf'OO 
Pl^H(0»-(tso<nr^r>.>otA^vOf^oooof^-,'^-0®orr(f^r-(i-t<nr»'Oooo\rio\r^<nf-t«o 


CO  H  H  ^  O'  -» 
M  i-( 

I  I 


N<0'»'»MN^H>»0'»<Nn'OBN'OinNO'inN>9ooi*)rtr.rtinoMfl 

r^i  r*4  r*4  H  H  H 

lilt  II  III  I  II  I 


ON'enN«aoi4in^nHN>»nu9'onnN(Mn«<r'4eoo''HNHNNNHs<«a 

HM<tNMn<no)n4inoNionNnNO>Aieininino'ono'Hn'0-«NHMininH<9 

OSN(SNNOlrtHO'NNOnrtri-JH'OO'O''fiNO00H»>OlfltOO"ONNNNinO'«N 

H  H  H  H  H  H  H 

II  I  I  I  I  I  I  III  I  I  I  I  I  I  I  I  I  I  I 


K9N«HnNN4«nn-oN«Hi/\oog\4onO'9'aaiHno'oat<oninHNoioie 

(MHHNHioON«oo'onoo-9ifi'-;insvt^coeooinmHOONin'OMoooiniflNHrt 

<*(sOHtONH<oco<nn«N«Hsrd'«HNio«H'Oio<»eo'00'Ooa(>>inNNnNiMn 

P'4  r-4  t-M  rH  *H  f-1  r— I  H  H  r- < 

II  I  I  I  I  I  I  II  I  I  I  I  I  II  I 


«n 

•  <C  H  |  rNPIHHNO'flO'HHNOOnOO'inO'O'IN'OH'O'OHHin'OOONOOHO'O 

O  4J  0<N<ns<PIP(nnNHNNnNWHMnNNNN(ONinNNr>NNNNNNlnnNN 
(031  OOOsrOGOOQOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 


r~fno'0'»OOOl*iP'OPr>'OP,,!srOOOOr^Orr>oorOf»imcoPir»Or'i^-9-OfnOP' 

09  N'O#Hpip!d«0»indN''nHpidN'O<,?I^NHO'pio'N8''OI,?»'»'<1''pJN»dH 
CM  r-4  H  M  H  N  NHNNHNHMrIHNNHH  rl  N  PI  rl  N  N  N  H  N  rt  N  H  H  N 

outoooooooom<-ioinotnooooo9oo>nooo90-<roooooopMOo 

3  H'0'CHiANnpiinNt09'<rov^<ONinrsaa'HNN«P4Nnntoo'<>i'»inHOioo 
Ninm'OiO'OininAinifiin'0«'einininininin'OiniOi'N'0'oeiO'ON'o>OP>ioiOP< 


m<6ioeininmioinin'SininNP>to«eiO'Ots'0>op>ioviN 


eneNO4eci4e'9PiH4NnNPion4ooaun9iOHNv9^innN0'NHo 

'9f'OM'IONO'0'Off'0<rHinPlHIPj'OP>0>0'ClH'PHOONeO®OMnONPIv': 

ONMOO»'HBH'»l,in  -  ‘  -  ‘ - ~ - - - -  - 


^HmnHni'orsO'O'ciHiPHiDNiaooQrmoNPi.; 

■»inOP''J'J-9NM,cic4'0'0(P'Or'H(NKi'ON<r>t(NO® 


'OO'POOOM'ONOM'O'IQOOOO'OO'OO-J'OOONMMNOCINON'J'OO'tO'TO 

r«>ervp4P«p.N'Of"orsrMB'Or>®r.NoO'0'or>p->ocoo'»oooioi''r"CMOP'ino' 


OQ'AMO'Of'BM't'O-jqoOOO-tOO-IM'J'P'O'IOO'tNOOMOOOOOnoN-J'OtM 

P'NN.p-POf'OPl'PHHrHnnNPIHNHnPINPlfPHNI'lN'jnrtH'inHCI 


«NSO'OHPin«in'0NB9iOHNn4in'eNe9O'QHPin<tin<ei>ao« 

HHHHHrlHHHHNNNMNNNNNNmW«rHP(n(n<n<*1(n 


152  88  .741  87.0  34.3  .029 


35 


The  parameter  estimates  obtained  via  the  two  methods  are  not 
strikingly  different.  The  only  obvious  difference  is  in  the  estimate 
for  the  fraction  of  life  term  for  both  systolic  and  diastolic  blood 
pressure  equations.  Table  6  also  includes  the  weight  assigned  to 
each  observation  by  the  SC  procedure.  Observations  4  and  1  received 
the  lowest  final  weights,  w^.  The  22-element  of  the  covariance 
matrix  is  quite  sensitive  to  changes  in  the  value  of  c  ant  the  resid¬ 
uals  of  Table  6  confirm  this. 

The  negative  coefficients  for  the  fraction  of  life  in  the  new 
environment  support  the  results  of  the  African  studies,  showing  that 
blood  pressure  decreases  with  the  length  of  time  in  the  new  environ¬ 
ment.  However,  it  appears  that  there  is  less  of  a  decrease  in  the 
diastolic  measurement  than  in  the  systolic  measurement.  As  expected, 
weight  has  the  effect  of  increasing  the  blood  pressure  measures,  but 
it  appears  to  have  a  greater  effect  on  the  systolic  measurement. 

The  four  term  linear  model  is  not  entirely  appropriate  since  there 
are  still  patterns  in  both  the  c  *  0  and  c  =  t  residuals.  In  this  case 
examination  of  the  c  =  0  residuals  would  have  led  to  the  two  potential 
problem  points,  1  and  4,  that  the  final  weights  w.  have  flagged 

lr  • 

for  further  scrutiny,  but  this  will  frequently  not  be  the  case. 

Our  procedure  may  be  directly  applied  to  multivariate  analysis  of 
variance.  A  generalized  inverse  may  be  used  in  place  of  the  inverse  of 
(7.5)  but  a  different,  more  direct,  approach  seems  preferable. 


ji1! 


ie 


8 .  Discussion 


He  have  presented  a  procedure  for  the  analysis  of  multivariate  normal 
data  which  is  at  once  very  general  and  easy  to  use.  He  have  successfully 
used  it  in  a  variety  of  applications.  It  has  proved  to  be  helpful  in  the 
identification  of  potential  outliers  and  in  the  process  of  model  evolution. 
The  procedure  will  tolerate  large  amounts  (more  than  100%)  of  contasdnation, 
provided  a  significant  portion  of  it  is  not  all  tightly  concentrated  in 
a  single  area  of  p-dimensional  space.  This  is  of  little  consequence  in 
any  event  because  the  final  weights  will  show  a  dramatic  sensitivity  to 
changes  in  the  index  c  and  the  source  of  the  sensitivity  should  always 
be  sought.  The  extent  to  which  the  index  c  may  differ  from  zero  is  a 
function  of  the  dimension  p  and  the  sample  size  n.  For  example,  if 
p  *  10  and  n  ■>  200,  the  procedure  is  likely  to  break  down  with  c  >  .15. 

It  is  highly  recommended  that  the  procedure  proposed  be  used  in  an 
exploratory  fashion. 
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